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ATMOSPHERIC CIRCULATION OF HOT JUPITERS: THREE-DIMENSIONAL CIRCULATION MODELS OF 
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ABSTRACT 

We present global, three-dimensional numerical simulations of the atmospheric circulation on HD209458b 
and HD 189733b and calculate the infrared spectra and light curves predicted by these simulations, which we 
compare with available observations. Radiative heating/cooling is parameterized with a simplified Newtonian 
relaxation scheme. Our simulations develop day-night temperature contrasts that vary strongly with pressure. 
At low pressure (< lOmbar), air flows from the substellar point toward the antistellar point, both along the 
equator and over the poles. At deeper levels, the flow develops an eastward equatorial jet with speeds of 3- 
4kmsec _1 , with weaker westward flows at high latitudes. This basic flow pattern is robust to variations in 
model resolution, gravity, radiative time constant, and initial temperature structure. Nightside spectra show 
deep absorption bands of H2O, CO, and/or CH4, whereas on the dayside these absorption bands flatten out or 
even flip into emission. This results from the strong effect of dynamics on the vertical temperature-pressure 
structure; the temperature decreases strongly with altitude on the nightside but becomes almost isothermal on 
the dayside. In Spitzer bandpasses, our predicted planet-to-star flux ratios vary by a factor of ~ 2-10 with 
orbital phase, depending on the wavelength and chemistry. For HD 189733b, where a detailed 8-/im light 
curve has been obtained, we correctly produce the observed phase offset of the flux maximum, but we do not 
explain the flux minimum and we overpredict the total flux variation. This discrepancy likely results from the 
simplifications inherent in the Newtonian relaxation scheme and provides motivation for incorporating realistic 
radiative transfer in future studies. 

Subject headings: planets and satellites: general, planets and satellites: individual: HD 209458b, methods: 
numerical, atmospheric effects 
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1. INTRODUCTION 

Of the ~ 255 extrasolar giant planets (EGPs) that have been 
discovered to date (see www.exoplanet.eu), the hot Jupiters 
have proved the easiest to characterize because of their high 
temperatures, sho rt orbital periods, and lik elihood of tran- 
siting their stars (Charbonne au et all 12007). The planetary 
mass and radius can be calculated from radial-velocity and 
transit data, allowing the surface gravity to be inferred; fur- 
ther more, the rapid spindown times for such close-in plan- 
ets dGuillot et alJl!99a) suggest that these planets are tidally 
locked, providing an estimate of their rotation rates. Such 
tidal locking implies that, in the absence of atmospheric 
winds, the dayside would be extremely hot and the nightside 
would be extremely cold. 

Impressively, infrared light curves have been obtained with 
the Spitzer Space Telescope for at least six hot Jupiters, 
helping to constrain the three-dimensional (3D) temperature 
pattern and atmospheric circulation of these objects. Con- 
tinuous 8-/j,m observations of HD 189733b over half an or- 
bit demonstrate that the planet's photospheric temperatures 
are relatively homogenized: the nigh tside is perhaps onl y 
~ 200-300 K colder than the dayside dKnutson et al.ll2007bl) . 
The exquisite temporal resolution allows detailed structure 
to be inferred. The light curve shows that the flux peaked 
16 ± 6° of orbital phase before secondary eclipse; a flux 
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map of the planet constructed from the light curve sug- 
gests that the hottest region lies ~ 30° longitude east of 
the substellar point and the coldest region lies a compara- 
ble amount west of the antistel lar point, providing cl ear evi- 
dence for advection by winds (Knutson et alJl2007bl) . 8-/im 
light curv es for HD209458b and 51 Peg b have also been 
obtained dCowan et al.ll2007h and similarly imply significant 
thermal homogenizationby winds, although these light curves 
contain insufficient temporal sampling to determine displace- 
ments (if any) of the hottest and coldest regions from the 
substellar and antistellar points. In contrast, Ups And b and 
HD 179949 exhibit large-amplitude phase variations, with no 
observationally significant phase shift, which suggests that 
these planets have large day-night temperature differences 
(perhaps > 500 K) that track the stellar in solation patterns 
dHarrington et alJl2006l: ICowan et afll2007l) . Together, these 
observations suggest that the degree to which winds homoge- 
nize the temperatures varies from planet to planet. 

There exist many additional observational constraints 
on the circulation, including dayside spectra, upper lim- 
its on albedo, and constraints on composition. For 
HD209458b, Spitzer secondary-eclipse photometry at wave- 
lengths of 3.6-24 /im sug gests the existence of a hot strato- 
sphere on the dayside dKnutson et alj l2007at iBurrows et alj 
l2007bl) . while HD 189733b seems to lack such a tem- 
perature inversion. Spitzer IRS spectra of these plan- 
ets' daysides from 7-14jnm lack obvious molecular ab- 
sorpti o n bands dGrillmair et al. 20071; iRichardson et alj|2003[ 
l2007t ISwain et aH |2007|). although deba te exists about 
the interpretation dFortnev & Marie vl 12007). Transit spec- 
troscopy observations ha v e yielded detections of sodium 
( Charb onneau et al. 2002; Redfiel d et alJ l2008h and water 
(Barman 2007; Tinetti et al. 2007) on both planets, and upper 
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limits on several other compounds. 

Several groups have investigated the atmospheric circu- 
lati on on hot Jupi t ers (f or a review o f the appro aches 
see IShowman et all [2007b ICho et all d2003l 120081) and 
Langton & Laughlin (2007, 2008) solved variants of the two- 
dimensional (2D) or quasi-2D equations on a rotating sphere. 
The advantage of this approach is that, by reducing the ver- 
tical resolution to one layer, simulations with high horizon- 
tal resolution can be performed. The idealizations implicit 
in such simplified models also allow the model dynamics 
to be relatively easily understood; a comparison of such re- 
sults with 3D models and observations can lend important 
insights into which dynamical processes cause which out - 
comes dVasavada & Showman 120051 IShowman et alj|2007l) . 
This has proved important, for example, in studying horizon- 
tal vortex/jet interactio ns on Solar-System giant planets and 
in Earth's stratosph ere dCho & Polvanill99fiiShowmanl200"7l: 
iPolvaniet al.ll 19951) . 

On the other hand, one-layer models do not represent the 
vertical structure of the flow, and they exclude inherently 3D 
processes such as baroclinic instabilities 6 , vortex tilting, and 
vertical wave propagation. In some contexts, such as the tro- 
pospheres of the Earth, Mars, Jupiter, and Saturn, 3D pro- 
cesses (e.g., baroclinic instabilities and convection) are cru- 
cial for determining the mean state because they convert po- 
tential energy to kinetic energy, hence causing the accelera- 
tions that pump the east-west jet streams. To some degree, 
such 3D processes can be parameterized in 2D models, but 
this introduces uncertainties regarding (for example) the rate 
at which the parameterized processes inject energy into the 
flow. In the hot- Jupiter context, a 3D approach is further moti- 
vated by the expectation that the radiativ e time constan t varies 
by orders of magnitude in the vertical dlro et al. 2005), which 
in the presence of the day-night heating gradient would cause 
patterns of temperature and wind that vary both vertically and 
horizontally in an inherently 3D manner. Because infrared 
spectra of planetary atmospheres depend sensitively on the 
vertical temperature profile, a full 3D representation of the 
winds and temperatures is also necessary for robustly predict- 
ing infrared spectra and wavelength-dependent light curves in 
Spitzer bandpasses. 

In light of the above, we have adopted a 3D approach 
to investigate the atmospheric circulation of hot Jupiter s 
dCooper & Showmanl2005L|2006tlShowman & Guilloti2 002). 
Coope r & Showman! (|2005[ 120061) and IShowman & Guillotl 
d2002l) focused on HD209458b and approximated the day- 
side heating and nightside cooling using a Newtonian cool- 
ing/heating scheme, which parameterizes the radiative heat- 
ing rate (in Ksec" 1 ) as (T ecj - T) / V ra( j, where r eq is the speci- 
fied radiative-equilibrium temperature profile (hot on the day- 
side, cold on the nightside), T is the actual temperature, 
and r m d is the radiative-equilibrium ti mescale, which was 
taken to be a function of pressure. In ICooper & Showman] 
d2005l). the ver tical structure of T eq and r ral j were taken from 
llro et al.l d2005). The simulations snowed development of sev- 
eral broad jets — including a superrotating 7 equatorial jet 

6 Baroclinic instabilities are a form of sloping convection that can occur in 
the presence of horizontal temperature contrasts in rotating, statically stable 
atmospheres. Lateral and vertical motion forces the cold air downward and 
the hot air upward, which releases potential energy that helps to drive the 
circulation. Most large-scale winter weather in the United States, Europe, 
and Asia is caused by baroclinic instabilities. 

7 SupeiTotation refers to eastward zonal winds, that is, winds moving faster 
than the planetary rotation. 



with speeds up to several km sec" — and exhibited horizon- 
tal (day -night) temperature di f ferenc es that varied strongly in 
height. Coope r~& Showmanl d2006l) extended these simula- 
tions to include carbon chemistry, which showed that inter- 
conversion between CO and CH4 should become chemically 
quenched at low pressure, leading to nearly constant abun- 
dances of these species everywhere above the photosphere 
at abundances that depend on the temperature in the deep 
(~ lObar) atmosphere. Using a multi-stream plane-parallel 
radiative-transfer code, iFortnev et alj d2006al) calculated the- 
oretical infrared light curves and spectra for HD209458b from 
the 3 D temperature patterns of Coo per & Showmanl (2005, 
2006). These calculations suggested that the day-night flux 
differences could reach factors of ~ 2-10 depending on wave- 
length and that the peak IR emission should lead the sec- 
ondary eclipse by ~ 2-3 hours in most Spitzer IRAC bands. 

Here, we present new global, 3D numerical simulations, 
spectra, and light curves that continue the rese arch program 
begun b y lCooper & Showmanl d2005Ll2006l) and lFortney et alJ 
(2006aJ). M ajor improvements ar e incorporate d in several 
areas. First, ]Cooper & Showmanl d2005l |2006) treated the 
day-night difference in T eq as a free parameter (with values 
ranging from 1 00-1000 K) rather than the result of a later- 
ally varying radiative-transfer calculation. This uncertainty in 
Ar eq directly translates into uncertainty in the light curves and 
spectra. Within the context of the Newtonian-cooling frame- 
work, the best approach — which we take here — is to ex- 
plicitly calculate the longitude, latitude, and height-dependent 
radiative-equilibrium temperatures for use in the Newtonian 
heating/cooling scheme. Ar eq is thus no longer a free param- 
eter but represents the true difference in radiative-equilibrium 
temperature from the dayside to the nightside. Likewise, we 
here calculate r rac j over the full 3D grid and express this as a 
function of both temperature and pressure rather than treating 
it as a function of p alone (bas ed on the ID global average 
tempe rature profile) as done by Cooper & Showman (2005, 

Second, we consider not only HD209458b (the canoni- 
cal planet emphasized by most previous circulation studies) 
but HD 189733b, and perform parameter variations in plane- 
tary rotation rate, gravity, and radiative time constant to de- 
termine how these parameters affect the behavior. Finally, 
we calculate infrared spectra and light curves from the 3D 
temperature patterns for comparison with available observa- 
tions, including the rec ent 8-/im light curve for HD 189733b 
dKnutson et alj l2007bl). secondary -eclipse photometry for 
HD 189733b dDeming et alJ [2006), and secondary-eclipse 
photometry for HD209458b, which now spans all availab le 
Spitzer channels dKnutson et alJl2007atrDeming et al.ll2005l) . 

A long-term goal is to couple the dynamics to a realistic 
representation of radiative transfer. However, this is a diffi- 
cult task, and experience with planetary general circulation 
models (GCMs) shows that the radiative-transfer calculation 
can become a computational bottleneck, making such a model 
computationally expensive. The computational and concep- 
tual simplicity of Newtonian cooling provide compelling ar- 
guments for documenting the extent to which a model driven 
by such simplified forcing can explain available observations. 
The results described herein thus provide a benchmark doc- 
umenting the best that can be done with this simplified ap- 
proach. 

In §2 we describe the dynamical model and present the cal- 
culations of laterally varying radiative-equilibrium tempera- 
ture and radiative time constant. §3 presents the basic dynam- 
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ical results and parameter variations. §4 describes the spectra 
and infrared light curves in Spitzer bandpasses as predicted 
by the simulations. §5 compares our results with pertinent as- 
pects of dynamical calculations by other authors. §6 shows 
that the global-scale flow of hot Jupiters should be close to 
local hydrostatic balance, and §7 concludes. 

2. MODEL 

2.1. Dynamics 

We solve the global, three-dimensional primitive equations 
in sph erical geometry using the ARIES/GEOS dynamical 
core (ISuarez & Takacslll995l) . The primitive equations are 
the standard equations for large-scale flow in stably strati- 
fied atmospheres whose horizontal dimensions greatly exceed 
the vertical dimensions. This is expected to be true on hot 
Jupiters, which have horizontal scales of 10 7 -10 8 mbut atmo- 
spheric scale heights of only 200-500km, leading to a hor- 
izontal:vertical aspect ratio of ~ 20-500. This large aspect 
ratio allows the vertical momentum equation to be replaced 
with local hydrostatic balance, meaning that the local vertical 
pressure gradient dp/dz is bal anced by the l ocal fluid weight 
pg (for a derivation see, e.g., Holton 2004). The primitive 
equations admit the full range of balanced motions, buoy- 
ancy (gravity) waves, rotationally modified (e.g., Kelvin and 
Rossby) waves, and horizontally propagating sound waves, 
but they filter vertically propagating sound waves. The hori- 
zontal momentum, vertical momentum, mass continui ty, and 
therm odynamic energy equations are as follows (e.g. jKalnavl 
UK)! p. 60-67): 



dx 

— = -V<J>-fkxv 
dt 

<9$ _ 1 
dp p 

_ duj 
V-v+— = 
op 



dT 
dt 



q lu 
c p pc p 



(D 
(2) 

(3) 

(4) 



where v is the horizontal velocity on constant-pressure sur- 
faces, u = dp/dt is the vertical velocity in pressure coor- 
dinates, $ is the gravitational potential on constant-pressure 
surfaces, / = 2£!sin0 is the Coriolis parameter, £1 is the 
planetary rotation rate (2ir over the rotation period), k is the 
local vertical unit vector, q is the thermodynamic heating 
rate (Wkg -1 ), and T, p, and c p are the temperature, den- 
sity, and specific heat at constant pressure. V is the hori- 
zontal gradient evaluated on constant-pressure surfaces, and 
d/dt = d/dt + x ■ V + uid/dp is the material derivative. Cur- 
vature terms are included in v • Vv. The dependent variables 
v, uj, p, and T are functions of longitude A, latitude <j>, 
pressure p, and time t . 

Note that the vertical velocity is nonzero — it enters via 
the mass-continuity equation and the thermodynamic energy 
equation. For a known atmospheric state at a given timestep, 
Eqs.[T]and|4]are integrated forward, leading to expressions for 
v over the 3D grid at the subsequent timestep. The horizontal 
divergence of these velocities is generally nonzero, and the 
vertical velocity u> is then evaluated for that timestep via Eq.[3] 
with use of the boundary conditions. For conditions relevant 
to hot Jupiters, the characteristic vertical velocity is typically 



~ 10msec" 1 near the photo sphere ([Showman & Guillotl2002l : 
ICooper & Showmanll2006L and see §6). 

The ARIES/GEOS dynamical core discretizes the equa- 
tions in longitude and l atitude using an Arakawa C grid 
(lArakawa & Lambl 1 19771) and adopts a pressure coordinate 
in the vertical. To maintain numerical stability, we follow 
standard practice an d apply Shapiro and polar filtering t o 
the time tendencies (Shapiro! 1 1970t ISuarez &Takacsi ri995). 
The top boundary condition is constant pressure and the bot- 
tom boundary condition is an impermeable surface, which we 
place far below the region of interest. These boundaries are 
free-slip in horizontal velocity. We solve the equations us- 
ing horizontal resolutions in longitude and latitude of 72 x 45 
or 144 x 90 with 30-40 layers spaced evenly in log p. The 
ARIES/GEOS model has been widely used in Earth and Mars 
studies and has been successfully b enchmarked against stan- 
dard test cases dHeld & S uarez 1993). 

For HD209458b, we adopt gravity, planetary radius, and 
rotation rate f2 of 9.81msec" 2 , 9.44 x 10 7 m, and 2.06 x 
10~ 5 sec _1 (implying a rotation period of 3.5 days), re- 
spectively. The corresponding values for HD 189733b are 
22.62msec" 2 , 8.22 x 10 7 m, and 3.29 x lO^sec" 1 (imply- 
ing a rotation period of 2.2 days), respectively. For both 
planets, c p = 1.23 x 10 4 Jkg _1 K _1 and the molar mass is 

2.36 x 10~ 3 kgmor', implying a specific gas constant of 
3523 Jkg" 1 K" 1 . We adopt the ideal gas equation of state. 

The top layer is placed at 0.7mbar. Radiative calcula- 
tions that match the planetary radius suggest that the tem- 
perature profiles converge to the interior adiabat at ~ 50 
bars on HD209458b and - 500bars on HD189733b (Fig.Q}. 
Accordingly, we place the bottom pressure at lOObars for 
HD209458b and 900bars for HD189733b. The intial condi- 
tions contain no winds and adopt a temperature profile result- 
ing from a ID global-average radiative-transfer calculation of 
the planet. Our basic results are insensitive to the initial con- 
dition, as is discussed more fully in §3. 

The heating rate is 



q T eq (\(j),p)-T(\,(j),p,t) 



(5) 



Cp Trad(/7, T) 

We assume that the obliquity and orbital eccentricity are zero, 
and that the planet rotates synchronously, which implies that 
the radiative-equilibrium temperatures T eq (X,(j>,p) and radia- 
tive time constants r ra d(/?,r) are time-independent; the fol- 
lowing subsection describes how we calculate them. The sub- 
stellar latitude and longitude remain fixed at 0°, 0° throughout 
the simulation. 

2.2. Radiative-equilibrium temperature structures 

We specify T eq in Eq. [5] by calculating the radiative- 
equilibrium temperatures as a function of longitude, lati- 
tude, and pressure for conditions relevant to HD209458b 
and HD 189733b. To do so, we calculate atmospheric 
pressure-temperature (p-T) profiles using the plane-parallel 
multi-stream r adiativ e transfer code first used for Titan by 
iMcKav et alj d!989l) and later extended to giant planets 
and brown dwarfs by Marley, Fortney, and collaborators 
(iMarlev et afl 11996] iBurrows et alJ [l997tjMarlev & McKavl 
Il999t iMarley etalJll999t iFortney et alJ 2005, 2006b aj). Al- 
though the metallicities of HD209458b and HD189733b are 
unknown, the host stars have near-solar metallicity. Here we 
adopt solar metal licity for the planets and assume local chem- 
ical equilibrium dLodders & Feglevll2002l 120061) . The calcu- 
lations use a large and continually updated opacity database 
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described in Freed man et al.l (120071) . For these exploratory 
calculations we neglect cloud opacity, although the "rainout" 
of elements that condense into clouds is always properly ac- 
counted for. 

For the present calculations, the opacities of TiO and VO 
are excluded, as has been standard for models of multi-Gyr- 
old hot Jupiters at the incident flux levels expected for gas 
giants beyond ~ 0.03 AU (for a Sun-like primary). One- 
dimensional atmosphere models suggest that gasous TiO and 
VO are cold trapped at pressures of ~ 10- lOObars where 
solid clouds incorporating Ti and V form, which suggests 
that TiO an d VO should be absent in the observable atmo- 
sphere (e.g. JHubeny et al.ll2003HFortney et al.ll2006bl) . How- 
ever, the rec ent detection of a ho t dayside stratosphere on 
HD209458b dKnutson et alj|2007al) suggests t hat future stud- 
ies of this planet should in clude TiO and VO (B urrows et alj 
l2007bHFortnev et al.l2007b . The limited da ta for HP 189733b 
are c onsistent with no TiO/VO opacity dFortnev & Marlevl 
120071) . 

The code treats thermal radiation from the planet from 
0.26-325 pm and incident radiation from the parent star from 
0.26-6.0yum. The calculation is performed at 60 model lay- 
ers from 0.3 mbar to 1000 bars, although the opacities at 
p > lOObars are quite uncertain owing to imperfect knowl- 
edge of the line broadening under these conditions. The cal- 
culation includes incident flux from the parent star, thermal 
flux from the planet's atmosphere, and thermal flux from the 
planet's interior, parameterized by Tjnt, the intrinsic effective 
temperature. 

To calculate radiative -eq uilibrium p-T profile s, we take an 
approach similar to that of iBarman et al.l d2005l) . At a given 
location on the planet's day side, the incident stellar flux ar- 
rives with an angle 8 from the local vertical. The cosine of 
this angle, /i, varies from 1 at the substellar point to at the 
terminator. Large p allows deeper penetration of stellar flux, 
greater absorption, and a warmer atmosphere. Smaller p gives 
a longer path length and airmass to a given pressure, implying 
shallower penetration of flux, less absorption (greater scat- 
tering), and a cooler atmosphere. We calculate p-T profiles 
for concentric rings of atmosphere that are symmetric around 
the normal at the subsolar point. The profiles are calculated 
at p = 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 
0.02, 0.01, and 0.0. The latter case, which is the radiative- 
equilibrium profile for an isolated object with the planet's 
known radius, gives the radiative-equilibrium profile across 
the entire nightside. All cases use a constant specific entropy 
at the bottom boundary, corresponding to the value needed 
to reproduce the planetary radius with an assumed heavy el- 
ement abundance of ^30-40 Eart h masses in the planeta ry 
interior, similar to that of Jupiter (Saum on & Guil lot 2004). 

Figs.QJ and b depicts the resulting profiles for HD209458b 
and HD 189733b, respectively. These atmospheric p-T pro- 
files show the atmosphere structure that each planet would 
attain under radiative equilibrium. The large temperature gra- 
dients implied on the planetary day side shows that vigor- 
ous dynamics is expected to occur. In chemical equilibrium, 
CO would dominate over CH4 across much of the dayside, 
but CH4 would dominate on the nightside. In the absence 
of any energy redistribution, the night sides of these plan- 
ets would be quite cold, with 100-mbar temperatures of 500 
K and 200 K for HD209458b and HD 189733b, respectively. 
The night sides abundances would be dominated by H2O, 
CH4, and NH3, and would be devoid of atomic alkalis, which 
would condense into clouds. Perhaps H2O cloud condensa- 



tion would occur as well. Dynamical redistribution of energy 
will of course alter the atmospheric temperature structure and 
chemistry. 

2.3. Radiative timescales 

We developed a simple method for calculating the temper- 
ature and pressure dependent radiative timescales across the 
large p-T space accessed by these atmospheres. We first cal- 
culate a suite of converged radiative-convective model atmo- 
spheres similar those shown in Fig.Q] Since we are interested 
in radiative time constants as a function p and T in these atm- 
sopheres, which may in principal be found anywhere in the 
atmosphere, at any value of fi, we calculated ~15 profiles 
across the p-T space shown in these figures, but for each pro- 
files assumed day-side average conditions, meaning /i=0.5. 
We increased and decreased the incident stellar flux to reach 
temperatures from ^100-2500 K for both planets. 

For each one of these profiles, we add a thermal perturba- 
tion 6T=10 K to a given model layer with pressure po- Be- 
cause this perturbed profile deviates from radiative equilib- 
rium, the perturbed layer emits a different flux than it receives, 
implying a cooling/heating that should cause the thermal per- 
turbation to decay. Given the perturbed profile, we perform 
a radiative-transfer calculation to calculate the net flux versus 
height. For any function fit), a characteristic timescale for 
variation off is f{df/dt)~ l (for exponential decay this would 
simply give the e-folding timescale). Thus, the characteristic 
timescale for decay of the thermal perturbation is ST divided 
by the heating rate expressed in Ksec" 1 , namely 



Trad = ST 



PCp 

dF/dz 



(6) 



where F is the net vertical radiative flux, p is density, z is 
height, and dF/dz is the volumetric heating rate ( Wm" 3 ). By 
varying po, we can calculate the r ral j as a function of pressure 
for every p-T profile in our suite of models. Tests with dif- 
ferent values of ST and those utilizing perturbations that span 
multiple model layers yield similar results. 

Fig. |2] shows the resulting r ra d for HD209458b and 
HD 189733b. These are surfaces on a regular p-T grid, in- 
terpolated from the original calculations. Some modest ex- 
trapolation was done on the high 771ow p and low 77high p 
corners. Perhaps the most striking aspect of the calculations 
is the strong pressure dependence, with r ra( j values changing 
by 8 orders of magnitude at a give temperature, from mil- 
libars to hundreds of bars in pressure. Temperature effects 
are also substantial. At millibar pressures we find a six or- 
der of magnitude change in r ra d across the temperatures that 
we explore, with the most significant temperature dependence 
coming below ^500 K, where r rac j becomes quite long. The 
calculations for HD 189733b extend to higher pressures than 
HD 209458b because of the expected deeper radiative zone 
in the atmosphere of HD 189733b. On average, the values 
for HD 209458b slightly exceed those for HD 189733b (by 
a factor of 1.5-3). This results from the lower gravity for 
HD 209458b, which leads to a "puffed up" atmosphere with 
greater mass per area across a given pressure interval. 

To our knowledge the only previous calculation of radiative 
time constants over a wide range in atmospheric pressure for 
hot Jupiters was bv llroetalj {2005) for HD 209458b. These 
calculations were limited to r rac j as a function of p only along 
the ID atmosph eric p-T profile that these authors calculated 
for the planet. Ilro et al.l d2005l) determined r ra d by adding 
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Gaussian temperature perturbations to their previously con- 
verged radiative equilibrium profile. They then performed 
a time-dependent radiative-transfer calculation to explicitly 
simulate the characteristic time needed for their model atmo- 
sphere to relax back to the radiative equilibrium structure. Our 
p-T-T m d surfaces shown in Fig. |2] cover a much larger phase 
space than was explored by these authors. Although our cho- 
sen methods are quite different, Fig. [3] s hows that o ur derived 
r ra< i values for the ID profile from llro et al.l (12005b are actu- 
ally very similar. Tabulated values of r ra d as a function of p 
and T for planets HD 209458b and HD 189733b are found in 
Tables 1 and 2. Nevertheless, the precise T ra d values should be 
viewed as uncertain because deviations from radiative equi- 
librium in the real atmosphere will not follow the idealized 
shape adopted here. 

3. RESULTS: DYNAMICAL SIMULATIONS 
3.1. Flow regime and dependence on parameters 

Co nsistent with the results of Coop er & Showman] ( 2005, 
2006), the imposed day-night heating contrasts lead to rapid 
development of large day-night temperature differences and 
winds reaching several km sec" 1 . Figures 14^5] show the tem- 
perature (greyscale) and horizontal winds (arrows) for our 
nominal simulations of HD 189733b and HD209458b at times 
of 900 and 800 days, respectively, by which time the temper- 
ature and wind patterns have reached a quasi-steady state at 
p < 1 bar. Panels are shown at pressures of 10 mbar, 100 mbar, 
and 1 bar, which bracket the range of pressures expected to be 
important for observable spectra and light curves. Both simu- 
lations use a resolution of 144 x 90 with 30 layers. 

As expected, the flow develops patterns of temperature 
and wind that vary horizontally and vertically in an inher- 
ently 3D manner. At low pressure (p < 10 mbar), the ra- 
diative time constants are shorter than the typical time for 
wind to advect across a hemisphere, and so the temperature 
patterns track the stellar heating — hot on the dayside and 
cold on the nightside, with approximate boundaries at the ter- 
minators (longitudes ±90°). The predominant flow pattern 
moves air away from the substellar point toward the antis- 
tellar point, both along the equator and over the poles. At 
deeper levels, the flow forms a banded structure dominated 
by an eastward equatorial jet extending from latitude ~ 30°N 
to 30°S (Figs. @H5]3 and c). In the range - 50-300mbar, 
the advection times are similar to the radiative times, lead- 
ing to large longitudinal temperature differences whose pat- 
terns are distorted from the day-night heating patterns; at 
p > 1 bar, however, the radiative timescale is longer than the 
longitudinal advection time and temperatures become homog- 
enized in longitude. The advection times in latitude exceed 
those in longitude, so nonzero latitudinal temperature dif- 
ferences persist even at relatively deep levels (Fig. HJ: and 
[5};). This banded structure results directly from the effects 
of planetary rotation, as has been shown by previous authors 
dShowman et al.l 120071; IShowman & Guiliol EOOl ICho et al l 
2003L 120081; IMenou et al.1 120031: ICooper & Showman! 120051: 
Dob bs-Dixon & Linll2008l) . 

Consistent with the r esults of previous studies 
(IShowman & Guillotl 12001 IMenou et all l200l IChoetal. 
20031 120081; lLangton & Laughhnl 120071: iDobbs-Dixon & Lin 
2008), the horizontal lengthscales of the simulated circulation 
patterns are comparable to a planetary radius. This results 
from the fact that the Rhines length and the Rossby deforma- 
tion radius, which tend to control the horizontal dimensions 
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FIG. 1. — Pressure-temperature ip-T) profiles for planet HD 209458b 
(top) and HD 189733b (bottom) as a function of fi. The hottest profile is 
li=1.0, and the profiles decrease in steps of 0.1 to /j=0,l. Profiles for /i=0.05, 
0.02, 0.01 are also shown, as well as the non-irradiated night side profile 
(p = 0). Condensation curves are shown as dotted lines and are labeled. Equal 
abundance curves for CO/CH4 and N2/NH3 are dashed. 

of the dominant flow structures, are comparable to a planetary 
radius for the modest rotation rates, hi gh temperatures, and 
large static stabilities relevant here (Showman & Guillotl 
2002; [Menou et al.ll2003t IShowman et alj|2007l) . 

The qualitative patterns of wind and temperature are similar 
in our HD 189733b and HD209458b simulations, suggesting 
that these planets may have comparable circulation regimes 
despite modest (factor of 1 .5-2) differences in gravity, rota- 
tion rate, and stellar flux. However, examination of Figs. |4] 
and [5] shows that the absolute temperatures differ substan- 
tially between these cases. Our HD209458b simulation pro- 
duce minimum and maximum temperatures of ~ 400 K and 
~ 1700-1800 K, respectively, in the observationally relevant 
layer from 10-1000mbar. In contrast, our HD189733b simu- 
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FIG. 2. — Radiative time constant (T ra( j) as function of pressure and tem- 
perature in the atmosphere of HD 209458b (top) and HD189733b (bottom). 
Time constants were derived from planet-wide average p-T profiles computed 
with the incident stellar flux increased and decreased to map out high and low 
temperatures. Small-scale bumpy structure in the surfaces is not physically 
significant. (See text.) 



lation is cooler, with temperatures ranging from 300 K to only 
~ 1400 Kin this same layer. Although HD189733b lies closer 
to its parent star (0.0313 AU as compared to 0.046 AU for 
HD209458b), the stellar luminosity is lower, explaining the 
cooler planetary temperatures. 

A careful comparison of Fig. Q] with Figs. |4]-|5] shows how 
the actual temperatures compare to the radiative equilibrium 
temperatures. For HD189733b, at low pressure (~ lOmbar), 
the local radiative-equilibrium temperature T eq ranges from 
155-1340K while the actual temperature ranges from ~ 300- 
1300K (fluctuating slightly in time). For HD209458b, at this 
same pressure, r eq ranges from 290-1670K while the ac- 
tual temperature ranges from 423-1687K. Both these com- 
parisons show that T lies very close to r eq on the dayside, 
but that T is substantially warmer than T eq on the nightside. 
This effect results directly from the temperature dependence 
of r ra d. At deeper levels, where r ra d is greater, T is substan- 
tially warmer than r eq on the nightside and cooler than T eq on 
the dayside. 

Figure [6] shows temperature profiles from our nominal 
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FIG. 3. — A comparison of the radiative time constant (r ra d) as function 
of pre ssure for the atmospheric p-T profile of HD 209458b fro m llro et all 
(2005). The dotted curve is the calculation from Iro et ID 120051) . The solid 
curve is an interpolation (using this same profile) in the p-T-T rl & surface for 
HD 209458b shown in Fig.0 
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FIG. 4. — Temperature (greyscale) and winds (arrows) for nominal 
HD 189733b simulation at 900 days. Resolution is 144 X 90. For this and 
all simulations in this paper, the substellar point is at longitude, latitude (0° , 
0°); the terminators are at longitudes ±90°. 
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FIG. 5. — Temperature (greyscale) and winds (arrows) for HD209458b 
simulation at 800 days. Resolution is 144 X 90. 
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FIG. 7. — Temperature (grayscale) and winds (arrows) for HD209458b 
simulation with twice the nominal Q at 800 days. Resolution is 144 X 90. 
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FIG. 6. — Temperature profiles for nominal HD209458b case averaged 
across a disk 45° in radius centered at the substellar point (solid), 90° longi- 
tude (dashed), antistellar point (dashed-dotted), and -90° longitude (dotted). 



HD209458b case. To better represent the broad-scale con- 
ditions relevant for spectra and lightcurves (to be discussed in 
§4), the profiles have been averaged across a circular region 
45° in radius centered at the substellar point (solid), longitude 
90° (dashed), antistellar point (dashed-dotted), and longitude 
-90° (dotted), all at the equator. The profiles demonstrate 
that dynamics exercises a zeroth-order influence on the verti- 
cal temperature structure. On the nightside (dashed-dotted), 
the intense radiative cooling leads to a temperature that de- 
creases strongly with height. As the equatorial jet transports 
this air eastward to the terminator at longitude -90° (dotted), 
deep regions (0. 1 to several bars) have continued to cool. In- 
tense stellar heating initiates as this air reaches the dayside, 
but this has an immediate effect only at low pressures where 
radiative time constants are short. A thermal inversion thus 
develops at low pressures (< 50mbar). Westward flow at 
high latitudes and at pressures < 10 mbar also contributes to 
this inversion, as such air comes from the hot dayside. By 
the time the superrotating equatorial jet transports the equa- 
torial air eastward to the substellar region (solid), heating has 
had a substantial effect down to pressures of several hundred 
mbar, leading to a deep thermal inversion layer. As the air 
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FIG. 8. — Temperature (greyscale) and winds (arrows) for HD209458b 
simulation with half the nominal f2 at 800 days. Resolution is 72 X 45. 



reaches the terminator at 90° longitude, rapid cooling aloft 
leads to low temperatures at pressures < 100 mbar, but the 
layer between 0.1-1 bar has longer radiative time constants 
and, having just crossed the entire dayside hemisphere, re- 
mains warm. The key point is that the T(p) structure varies 
strongly across the globe and deviates strongly from the pre- 
dictions o f ID radiative- e quilib rium models, as previously de- 
scribed in iFortnev et al.l d2006af) . The day-night temperature 
differences are ~ 300 K at lbar, ~ 600 K at 100 mbar, and in- 
crease to ~ 1000K at the top of the model. 

To test the sensitivity of our results to model resolution, we 
also performed simulations of HD209458b at a lower horizon- 
tal resolution of 72 x 45 with 40 layers. These simulations 
produced patterns of wind and temperature that are nearly 
identical to Fig. [5] While we cannot rule out that the behav- 
ior could change at a resolution higher than considered here, 
this test suggests that, within the context of our forcing ap- 
proach, our horizontal resolutions are sufficient to resolve the 
global-scale flow. This makes sense in light of the fact that 
(i) the predominant forcing consists of heating gradients with 
a large (hemispheric) length scale, and (ii) the Rhines length 
and Rossby deformation radius, which determine the predom- 



inant horizontal widths of jets and baroclinic eddies, are com- 
parable to the planetary radius. Our resolutions are sufficient 
to resolve these lengthscales. This differs from Jupiter, where 
a relatively fine horizontal resolution is required to resolve the 
small deformation radius of ~ 2000 km. 

To determine the influence of planetary rotation rate on the 
flow geometry, we performed HD209458b parameter varia- 
tions with double and half the nominal rotation rate (rotation 
periods of 1.75 and 7 days, respectively); the results are de- 
picted in Figs. [7] and [8] These models are otherwise essen- 
tially identical to the nominal HD209458b case, including the 
use of an insolation pattern relevant for a synchronously ro- 
tating planet (i.e., the substellar point is locked at 0° longi- 
tude). Although the basic flow regimes are similar to those 
obtained in the nominal case (Fig. |5]), a careful comparison 
reveals several key differences. Most importantly, as can best 
be seen at ~ lOOmbar, the superrotating equatorial jet is nar- 
rower in the high-rotation-rate case (Fig. |7]i and wider in the 
low-rotation-rate case (Fig. [8]) as compared to the nominal 
simulation. The midlatitude flow is, on balance, weaker at the 
higher rotation rates. Furthermore, the transition from sim- 
ple dayside-to-nightside flow to a banded jet pattern occurs at 
lower pressures at the higher rotation rates. Fig.|5h and|7^ de- 
velop equatorial flow that is eastward at most longitudes; on 
the other hand, Fig.[8h exhibits eastward flow at longitudes 
to 180° and westward flow at longitudes -100 to 0°, leading 
to a pattern nearly symmetric in longitude about the substel- 
lar and antistellar points. As a result, substantial differences 
arise in the temperature and wind patterns at p < 100 mbar in 
the nightside quadrant west of the antistellar point (longitudes 
90-180°). 

These differences in the jet patterns can be better appreci- 
ated by examining the plots of zonally averaged zonal wind 
shown in Fig. [9] The top panel gives the nominal HD209458b 
simulation (as in Fig. [5), and the middle and bottom pan- 
els depict the high- and low-rotation-rate parameter varia- 
tions (as in Figs. ITJjH). The equatorial jets in the nominal 
and low-rotation-rate cases have latitudinal widths (character- 
ized by full-width-at-half-max) approximately 20% and 40% 
greater, respectively, than in the high-rotation-rate parame- 
ter variation. This agrees qualitatively with the expectation 
that, for constant wind speed, scale height, and static stabil- 
ity, faster rotation implies smaller Rhines length and Rossby 
deformation radius — and thus should lead to narrower jets. 
On the other hand, the simulated variations in jet width are 
less than suggested by Rhines scaling, which predicts that, 
at constant wind speed, factor-of-four variations in rotation 
rate should produce factor-of-two variations in jet width. The 
fact that the forcing lengthscale is similar to the jet scale and 
that it remains constant throughout may play a role in muting 
the rotation-rate sensitivity. Interestingly, the jets penetrate 
slightly deeper in the higher-rotation-rate cases. The absence 
of the eastward jet at p < 10 mbar in the low-rotation-rate case 
(Fig. [9};) results from the fact that the flow is largely symmet- 
ric in longitude about the substellar point at these low pres- 
sures; the eastward and westward branches of the circulation 
cancel out in a zonal average, leading to weak zonal-mean 
wind speeds. 

We also performed HD209458b simulations with half and 
double the nominal gravity (5 and 20msec" 2 ). These simu- 
lations are very similar to the nominal case shown in Fig. |5j 
although the equatorial jet speed in the low-gravity case is 
slightly weaker than in the nominal and high-gravity cases. 

To determine the sensitivity to the radiative time constant, 
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FIG. 9. — Zonal-mean zonal velocity versus latitude and pressure for 
HD209458b simulations at 800 days with nominal rotation rate (top), double 
the nominal rotation (middle) and half the nominal rotation (bottom). Scale 
bar shows wind speed in msec -1 ; positive is eastward and negative is west- 
ward. The equatorial jet is wider when the rotation rate is lower. 

we performed an HD209458b simulation analogous to Fig. [5] 
but multiplying r ra d by a factor of 10 at all temperatures and 
pressures. This case is particularly relevant because of the un- 
certainty in our calculated T ra d- As expected, the winds in this 
simulation develop more slowly and the day-night tempera- 
ture differences are more muted than in the nominal cases. 
Qualitatively, however, the overall circulation pattern strongly 
resembles that in our nominal cases. 

We wish here to correct the recor d regarding the is- 
sue of sensitivity to initial conditions. lLangton & La ughlin 
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FIG. 10. — Heat source, expressed in Ksec~', associated with day-to-night 
horizontal advection across the terminators for our nominal HD189733b sim- 
ulation. Dashed and dotted lines show heat transport across the dawn and 
dusk terminators, respectively (third and fourth terms on the right side of 
Eq.[5J. Solid line shows their sum. 



d2008l) speculated that the sim ulated wind speeds reported 
in ICoope~& Showman (2005) occurred because their simu- 
lations were initialized with the nightside temperature rather 
than the global-average temperature. This is incorrect. Our 
simulation results are insensitive to the initial temperature 
profile, particularly at pressures < lObars; cold and hot ini- 
tial tempe ratures lead to very simil ar end states. Simula- 
tions by Coop er & Showman! d2006l) explicitly demonstrate 
this fact: they initialized their simulations with the global- 
average temperature profile rather than the nightside tem- 
perature, yeMhejr_simuJation result s are nearly identical to 
those in ICooper & Showmanl (120051) . This insensitivity oc- 
curs because the day-night temperature pattern that develops 
in the simulation (which via Eq. [redetermines the pressure- 
gradient forces that drive the winds) is dominated by the radia- 
tive heating/cooling patterns and hence quickly loses memory 
of its initial condition. 

3.2. Day-night heat redistribution 

Several recent ID radiative-transfer models include param- 
eteriza tions of the day-night heat r edistribution by the circu- 
lation (Bur rows et al.ll2007a[ 120061). In these studies, an ar- 
bitrary fraction of the absorbed starlight is removed from the 
dayside and added to the nightside. To guide such efforts, 
we here quantify the magnitude of redistribution that occurs 
in our simulations. The goal is to obtain an equation for the 
mean nightside temperature and its modification by day-night 
advection. Starting with the thermodynamic energy equation 
we express the advection terms in divergence form and 
invoke the continuity equation ([3]l to cancel out terms, yield- 
ing 

dT _ , d , q u 



a+v .(vr) + -( W r ): 



PCp 



(7) 



We then horizontally average the equation over the nightside 
hemisphere, yielding 

9 (T) night 

dt 



P I night 



+ ( — -7r^) "(V-(v7: 
\P c p dp / rag ht 



night 



(8) 
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where (X) n i g h t = (2ira)~ l f X dA, a is the planetary radius, the 
integral is taken over the nightside hemisphere (of area 27ra), 
and dA is the increment of horizontal area on the sphere. Us- 
ing the divergence theorem, we can then rewrite (V • (vr))ni g ht 
as the product of the zonal (east-west) wind and temperature 
integrated along the day-night terminators: 



HD 209458b, 800 days 



9 (T) night 



dt 



u> d 



+ / — -—(u;T) 

\P C P d P /night 
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2na 
1 

2na 



C P I night 
90° 

"( Adawn, ^, /?)r( Adawn , 0, /?) 

90° 
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(9) 



"(Adusk, Adu S k, <j>,p)d(j) 



Here, Adawn and Ad us k refer to the longitudes of the "dawn" 
and "dusk" terminators (-90° and 90°), and we have used 
the fact that for zero obliquity the terminators lie at constant 
longitude. 

Equation [9] provides a framework for quantifying the day- 
night heat redistribution. On the right side, the first and sec- 
ond terms in brackets are the temperature change due to ra- 
diative heating/cooling and vertical motion, respectively. The 
final two terms are the "heat source" associate d with thermal- 
energy tran sport from dayside to nightside. iBurrows et al.l 
(2007a, 2006) took the approach of introducing an ad hoc pa- 
rameterization for the third and fourth terms on the right side 
of Eq. [9] Here, we explicitly calculate these terms from our 
simulations. 

Figure [10] shows these heat transports versus pressure for 
our nominal HD 189733b simulation. The dashed and dotted 
curves show the transport across the terminators at longitude 
-90° and +90°, respectively, and the solid curve shows their 
sum. Several points deserve mention. First, the individual 
transports (dashed and dotted curves) are large at low pres- 
sure and decay toward zero at pressures > 100 bars, which 
results from the fact that, in these simulations, the circulation 
is strong aloft and dies out with depth. 

Second, these individual terms (dotted and dashed curves in 
Fig. [Tol l sum constructively at pressures < 50mbar but have 
opposite signs at greater pressures. This results from the fact 
that, at low pressure, the circulation transports air from day- 
side to nightside across both terminators (see Figs.|U|5]l, lead- 
ing to a positive heat source for the nightside across both ter- 
minators. In contrast, at depth, the superro taring equatorial 
jet transports air from the dayside to nightside across the dusk 
terminator (leading to a positive contribution in Fig.[l0] dot- 
ted) but transports air from the nightside to dayside across the 
dawn terminator (leading to a negative contribution in Fig.fTOl 
dashed). At the altitudes of the superrotating jet, the day-night 
heat transports across the two terminators largely cancel. As a 
result of these effects, the net advective heat source is almost 
zero at pressures exceeding ~ 1 bar but becomes large at lower 
pressure (solid curve). The magnitudes reach ~ 0.02 Ksec" 1 
near the top of the model. 

The magnitude of the net advective heat source (solid curve, 
Fig. [Tol l increases monot onically with a l titude. This differs 
from the scheme used by [Burrows et al.l (l2007bllah to predict 
light curves and secondary-eclipse depths, which confined the 
heat transport to pressures between ~ 0.01 and 0.1 bars. Our 
simulations suggest that a continuous, upward-increasing pro- 
file of temperature modification may be preferable. 

Interestingly, pressures where the advective temperature 
modification is greatest (i.e., the top of the model) have 
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FIG . 11. — Emergent flux density ( erg sec - ' cm -2 Hz~' ) from the nominal 
HD209458b simulation at six phases assuming equilibrium chemistry: black 
— nightside, as seen during transit; red — 60° after transit; green — 120° 
after transit; dark blue — dayside, as seen during secondary eclipse; light 
blue — 60° secondary eclipse; and magenta — 120° after secondary eclipse. 
The key in the upper right corner is color-coded with the spectra to illustrate 
the sequence. For reference, the dotted blue curve shows the spectrum of a 
ID radiative-equilibrium model atmosphere. Thin dotted black lines at the 
bottom of the figure show normalized Spitzer band passes and the letters at 
the top show locations of the H, K, L, and M bands. 
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FIG. 12. — Light curves versus orbital phase calculated in Spitzer band- 
passes for HD 189733b. From top to bottom, the light curves are for wave- 
lengths of 24 fim (magenta), 16/^m (light blue), 8fj,m (dark blue), 5.8 /im 
(green), 4.5 fim (red), and 3.6 fim (black), re s pective ly. Overplotted is the 
Spitzer 8-/tm light curve from Knutson et al. (2007b) in dark blue and the 
16-/im secondary-eclipse depth from Deming et al. (2006) in light blue. 



the largest day-night temperature difference. The enormous 
radiative heating/cooling rates at these levels overwhelm 
the transport, allowing near-radiative-equilibrium conditions 
(with large day-night temperature differences) to be main- 
tained. In contrast, at deeper levels, the transport modifies 
the temperature more weakly, but the heating/cooling rates 
are low, allowing relatively homogenized day-night tempera- 
tures to be maintained. This situation belies simple descrip- 
tions (common in the literature) of planets with large day- 
night temperature differences having "inefficient" redistribu- 
tion and those with homogenized conditions as having "effi- 
cient" redistribution. 

4. SPECTRA AND LIGHT CURVES 
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FIG. 13.— Same as in Fig,fl2lbut for HD 209458b. Diamonds give mea- 
sured Spitzer secondary-eclipse depths with 1-cr error bars. Wavelengths of 
measurements and light curves are color coded to each other an d to the num- 
bers given in the plot. Magenta is 24 fim iDeming et aL 12005) for the lower 
point and D. Deming, personal communication, for th e upper point), light 
blue is 16/^m, dark blue is 8/im iKnutson et al. 2007b), green is 5.8 fim, red 
is 4.5 fim, and black is 3.6 ^m (all from Knutson et al. 2007a). 
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FIG. 14. — Planet-to-star flux ratio versus wavelength for our nomi- 
nal HD209458b simulation at the time immediately before/a fter secondary 
eclipse. Points show mea sured secondary-eclipse depths I De ming et al.1 
I2005t IKnutson et"al]|2007bl lal). As discussed in the text, our inability to fit 
the points at 4.5 and 5.8 fim probably results from the lack of a stratosphere 
in our simulations. 

We calculate spectra and light curv es from our simulation s 
following the methods described in Fortney et al. (2006a). 
The 3D temperature field from a dynamical simulation at 
a given time can be viewed as numerous individual Tip) 
columns, one for each (longitude, latitude) location on the 
grid. Our low- and high-resolution simulations contain 3240 
and 12,960 such columns, respectively. At any orbital phase 
(except during secondary eclipse), half of these columns will 
be visible from Earth. Assuming solar metallicity, we run 
each such profile through our radiative-transfer solver (§2.2) 
to determine the net emergent flux density versus wavelength 
for that column. This calculation properly includes the appro- 
priate value of /i0 for each column, where is the cosine of 
the angle between local vertical of each column and the line 
of sight to Earth; this naturally accounts for any limb bright- 



ening or darkening. The emergent flux calculated for each 
column is then weighted by the apparent area of that patch of 
the planet as viewed from Earth; these are then summed to 
obtain the total emergent flux densi ty of the planet at a given 
orbital phase. iFortney et al.1 (l2006al) gives further details on 
the method. 

Fig. QT| depicts the resulting spectra for our nominal 
HD209458b case under the assumption that local chemical 
equilibrium holds; qualitatively similar behavior occurs for 
HD 189733b. The spectra vary strongly with orbital phase. 
During secondary eclipse (black solid curves), the emergent 
flux is low and deep absorption bands of CH4 and H2O ap- 
pear. Away from transit, as the dayside comes into view, the 
infrared flux increases and the absorption bands become shal- 
lower. Near secondary eclipse (dark blue), the dayside faces 
Earth and the spectra exhibit only modest deviations from a 
blackbody. Interestingly, at this phase, the weak spectral fea- 
tures that exist have flipped into emission. 

These phase variations in spectra result directly from the 
fact that the simulated Tip) structure depends strongly on lon- 
gitude. The radiative time constant is short at low pressure and 
long at high pressure (Pig. O. Thus, on the nightside, the air 
aloft cools more rapidly than the air at depth, leading to a tem- 
perature that decreases strongly with altitude (Fig. |6j. Con- 
versely, on the dayside, the air aloft warms more rapidly than 
the air at depth, leading to a shallower structure that is quasi- 
isothermal or even exhibits a temperature inversion where T 
increases with altitude (Fig. [6}. This latter behavior is the 
reason for the weak emission features in the dayside spectra 
(dark blue curves in Fi g . ITTb . 

We now t urn to infrared light curves. Integrating our spec- 
tra and the iKuruczl d 19931) model of the star HD189733 or 
HD209458 (as appropriate) over the Spitzer bandpasses, we 
calculate the planet-to-star flux ratios versus orbital phase. 
Figures [12] and Qj] display these light curves for equilib- 
rium chemistry for our nominal models o f HD 189733b and 
H209458b. As in IFortney et all (l2006al) . the light curves 
here show a phase offset that results from the distortion of 
the temperature field by dynamics — most importantly, the 
eastward displacement of the hottest region from the sub- 
stellar point and coldest region from the antistellar point, 
which causes the flux minima and maxima to occur be- 
fore transit and secondary eclipse, respectively. The off- 
set is greatest (~ 30°) at 3.6/zm and smallest (~ 20° for 
HD189733band 13° for HD209458b) at 24 fim. As described 
in lFortnev et al.l d2006al) . this wavelength-dependence results 
from the pressure-dependence of the photospheres, which 
tend to be at lower pressures for the longer Spitzer wave- 
lengths. In our simulations, lower pressures have temperature 
patterns that more closely track the stellar heating patterns 
(Figs.[U-[5]l, leading to the smaller phase shifts. 

In our simulations, the infrared light curves reach a steady 
pattern relatively quickly. Light curves calculated at 100 days 
are very similar to those calculated at 900 days. 

Fig-Q~2]compares our light curves to the 8-/xm Sp itzer IRAC 
light curve of HD189733b (IKnutson et al.ll2007bl) . Interest- 
ingly, our simulations correctly produce the observed off- 
set of the flux maximum from the time of secondary eclipse 
(16 ± 6°); however, we do not explain the location of the flux 
minimum, which occurs before transit in our simulations but 
after transit in the observations. To explain this feature would 
require a cold region to the west of the antistellar point. Be- 
cause the low-latitude winds in our simulations are predomi- 
nantly eastward, the cold region in our simulation is displaced 
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to the east rather than the west. In the context of the simula- 
tions presented here, simultaneously explaining both an east- 
ward displacement of the hot region and a westward displace- 
ment of the cold region is difficult. 

If the circulation is quasi-steady, one possibility is that the 
dayside and nightside 8-/im photospheres sample different 
pressures; the observed light curve might be explained if the 
equatorial jet speed were eastward at the pressure of the day- 
side photosphere and westward at the pressure of the nightside 
photosphere. Another possibility is that a localized region 
of strong ascending motion (perhaps driven by absorption or 
breaking of atmospheric waves) could transport low-entropy 
(hence low-temperature) air to the photosphere level in a lo- 
calized region, perhaps explaining the low-flux region. Al- 
ternatively, perhaps meandering hot and cold vortices could 
lead to time-var i able temperature patter ns and light curves 
dCho et alJ 12001 iRauscher et all l2007bllah : the low-flux re- 
gion might then result from fortuitous placement of a cold 
vortex to the west of the antistellar point at the time of the 
observations. However, we caution that the star HD 189733 
exhibits starspots that could influence the light curve; it is un- 
clear whether all the structure in the light curve results from 
the planet rather than variability in the star. Future observa- 
tional and theoretical work may resolve this ambiguity. 

Perhaps more importantly, our simulations greatly under- 
predict the 8 /im flux on the planet's nightside (Fig.fT2l). sug- 
gesting that our nightside temperatures are too cold at the 8- 
pm photosphere. At the expected 8-/im photospheric pressure 
of ~ 30mbar, our radiative time constants vary from ~ 10 4 sec 
at 1500K to 10 5 sec at 500 K and rapidly rise to 10 7 sec at 
temperatures less than 500K (Fig. |2). For an equatorial jet 
speed of ~ 3 km sec" 1 relevant for our simulations, the time 
to advect air parcels across a hemisphere is 10 5 sec; thus, in 
our simulations, the nightside temperature at the %-pm photo- 
sphere rapidly cools to ~ 500 K, at which point further cool- 
ing is inhibited by the rapidly increasing T ra d. In contrast, 
the nightside 8-/zm b rightness temperature inferred from the 
IKnutson et all d2007bl) light curve is 973 ± 33 K. 

The fact that we overpredict the 8-/im day-night flux vari- 
ation (and underpredict the absolute nightside flux) prob- 
ably results from the simplifications inherent to the New- 
tonian heating/cooling scheme. Newtonian heating/cooling 
neglects the fact that the radiative equilibrium temperature 
and timescale can depend on the atmosphere's dynamical re- 
sponse, particularly when actual temperatures are far from ra- 
diative equilibrium. In reality, nonlinearities not included in 
the Newtonian relaxation framework can lead to radiative in- 
teractions between vertical levels that are not accounted for 
here. For example, because of the slow radiative heating rates 
at deep levels, the temperatures at p > 1 bar become homog- 
enized between day and nightsides; upwelling infrared radia- 
tion from this level will then warm the entire p-T structure at 
p < 1 bar and lead to nightside radiative-equilibrium temper- 
atures that — because of this dynamical response — exceed 
those in Fig. Q] This could lead to warmer nightside condi- 
tions than predicted here. 

Next consider HD209458b (Fig. H"3"V From several brief 
Spitzer observations at different phases, ICowan et alJ (120071) 
place a 2-er upper limit of 0.0015 on the peak-to-peak flux 
variation at 8 ^m. In contrast, our light curve calculated un- 
der the assumption of equilibrium chemistry exhibits a peak- 
to-peak variation at 8 /im of ~ 0.0020, modestly higher than 
the upper limit. Regarding secondary-eclipse photometry, 
we match well the observed eclipse depths at 3.5 and 8/im 



( IKnutson et al.ll2007al) . Two meas urements exist of the 24- 
lim eclipse depth, 0.0026 ± 00046 (iDeming et alJl2005l) and 
0.0033 ± 0.0003 (D. Deming, personal communication); our 
simulated light curve is marginally consistent with the lat- 
ter. On the other hand, we underpredict the 4.5-pm and 5.8- 
pm eclipse depths by ~ 40% (see also Fig. fT4T> . The high 
flux at these wavelengths has been interpret ed as resulting 
from a thermal inversion at pr essures < 0. 1 bar (IKnutson et alJ 
l2007at IBurrows et al.ll2007bh . Producing this feature may re- 
quir e the existence of stratospheric absorbers su ch as TiO and 
VO dFortnev et all 120071; IBurrows etaill2007bllah . which are 
not included in the present models. 

Although in chemical equilibrium the primary carbon car- 
rier would be CO on the dayside and CH4 on the nightside, 
Cooper & Showman (2006) showed that interconversion be- 
tween CO and CH4 should become chemically quenched at 
low pressure, leading to nearly constant abundances of these 
species everywhere above the photosphere. For HD209458b, 
their simulations suggested a quenched CH4 mole fraction 
corresponding to ~ 1-2% of the total carbon. Although 
ICooper & Showmanl d2006l) did not consider HD 189733b 
specifically, their "cold" parameter variation has similar tem- 
peratures to our HD 189733b simulation; the quenced CH4 
abundance was ~ 20% of the total CO+CH4 in this case. The 
specific values should be viewed as quite uncertain. 

We thus recalculated the phase-dependent spectra and 
lightcurves assuming constant CO and CH4 mole fractions 
rather than chemical equilibrium. Having only a fraction of 
the carbon in CH4 on the nightside (compared to the equi- 
librium case where 100% of the nightside carbon resides in 
CO) lessens the depth of the nightside absorption features at 
3.6 and 8/im, making the nightside brighter at these wave- 
lengths. 8 The effect is modest when the fraction of carbon in 
methane is 20% but strong when it is only 2%. In the latter 
case, the day-night flux variation at 3.6/im is reduced from 
a factor of ~ 10 to ~ 2 w hile that at 8 pm is re duced from a 
factor of ~ 4 to ~ 2.5 (see lFortney et alJl2006al Fig. 9). 

By lessening the predicted day-night flux variation, this dis- 
equilibrium effect allows us to fit the upper limit of 0.0015 
on the p eak-to-peak planet /star flux variation of HD209458b 
at 8/im dCowan et alJl2007l) . It could also potentially help 
provide a better fit to th e S-p,m light curve of HD 189733b 
(Fig. Q21 IKnutson etail l2007bl) . but only if the planet's 
methane abundance is relatively small (perhaps 5% or less 
of the total carbon if the metallicity is solar). Clearly, light 
curves in the 3.6/im Spitzer band (and the 4.5/im/3.6/im 
band ratio, which is sensitive to the CH4/CO ratio) will pro- 
vide important constraints on possible disequilibrium chem- 
istry. 

5. COMPARISON TO OTHER WORK 

Showman et alJ d2007h reviewed the various approaches 
adopted for investigating the atmospheric circulation of hot 
Jupiters; here we expand up on key points. 

Of the published models, Dobbs- Dixon & Lir] (12008) most 
strongly resembles the approach taken here. They performed 
3D numerical simulations of the Navier-Stokes equations in 
the low- and mid-latitudes (excluding the poles). Radiation 
was parameterized by flux-limited diffusion using Rosseland- 
mean opacities. They parameterized stellar irradiation by im- 
posing a spatially varying temperature at the upper bound- 

8 In our simulations, the dayside fluxes are less sensitive to the CH4/CO 
ratio because of the near-isothermal dayside conditions. 
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ary (hot on the dayside and cold on the nightside); down- 
ward diffusion of energy at the model top thus induces heat- 
ing on the dayside while upward diffusion out the model top 
induces cooling on the nightside. Their nominal simulations 
adopt opacities relevant to interstellar-medium (ISM) gas and 
grains; as they point out, this is probably an overestimate for 
planetary atmospheres, where significant particle growth and 
settling is expected. They thus also performed cases with 
opacities reduced by factors of 100 and 1000 relative to the 
ISM and suggest that 100-fold reduction may be most appro- 
priate for hot Jupiters. 

For diffusion to be a good approximation, the mean-free 
path of the relevant photons must be much less than the 
length scales over which the atmospheric properties vary ver- 
tically. This is an excellent approximation in the deep interior 
where opacities are high. However, the diffusion approxima- 
tion breaks down near the photosphere, where a large frac- 
tion of infrared photons can escape directly to space. This 
is the important region for infrared spectra, light curves, and 
secondary-eclipse depths. In this region, radiative fluxes can- 
not be represented as a diffusivity times the local temperature 
gradient, as postulated in the diffusion approach; instead, the 
radiative transfer involves nonlocal inter actions between dif- 
ferent atmospheric levels. Nevertheless. iDobbs-Dixon & Linl 
(2008)'s approach, like ours, serves its intended purpose of 
providing an approximate, computationally efficient means to 
force a plausible flow. 

A further difference between IDobbs-Dixon & Linl ((2008) 
and the present study involves the representation of dynam- 
ics. They solved the Navier-Stokes equations rather than the 
primitive equations; we return to this issue in §6. Addition- 
ally, they included the centrifugal acceleration in the equation 
of motion; this differs from the approach taken in most plane- 
tary circulation models, which typically account for the gravi- 
tational relaxation of the planetary interior to the planetary ro- 
tation by absorbing the centrifugal force into the gravity. 9 At 
the present time, however, the greatest source of uncertainty 
lies in the treatment of radiation rather t han dynamics. 

Des pite the differences between Dob bs-Dixon & Linl 
(2008)'s approach and that adopted here, their basic results 
share strong similarities with ours. Because their case explor- 
ing opacities 100 tim es less than ISM values see ms most rele- 
vant to hot Jupiters (IDobbs-Dixon & Lirill2008l) . this is their 
best simulation to compare with ours. In agreement with 
our results, this case develops a banded pattern at the pho- 
tosphere with a broad, planet-encircling superrotating equa- 
torial jet reaching speeds of 3kmsec _1 and weaker (~ 1— 
2 km sec -1 ) westward jets in the midlatitudes (see top panels 
of their Figs. 9-10). Their day-night temperature differences 
at the photosphere reach ~ 500 K. Although our simulated 
day-night temperature difference at the photosphere differs 
somewhat from theirs (~ 800K and ~ 400K for our nominal 
and 10 x nominal r ra( j cases, respectively), the qualitative sim- 
ilarities between their simulation and ours are striking given 
the distinct approaches. 10 These similarities, which lend cre- 

9 For example, a hypothetical rotating planet with no radiative forcing, 
no initial horizontal temperature/pressure gradients (i.e., along the rotation- 
ally modified equipotentials), and zero initial winds should remain motion- 
le ss over time. If the initial c ondition is spherically symmetric, as assumed 
by Dobbs-Dixon & Lin 1 2008), then inclusion of the centrifugal acceleration, 
however, would cause equatorward motion in such a case — effectively forc- 
ing the equatorial bulge to develop in the modeled thin atmosphere rather 
than throughout the planetary interior. 

10 Their cases using ISM opacities exhibit some qualitative differences in 



dence to both approaches, are encouraging because they sug- 
gest that the basic flow regime on hot Jupiters is relatively 
in sensitive to differences in m odel formulation. 

Langton & Laughlir] d2007l) performed numerical simula- 
tions of the atmospheric flow on hot Jupiters using the one- 
layer shallow-water equations, which govern the behavior of 
a thin layer of hydrostatically balanced, constant-density fluid 
on a sphere. The m omentum and continuity equations are 
re.g. JPedloskvl[T987l Chapter 3) 

^-=-g\7h-fkxy (10) 
at 

^=-V-(/w) (11) 
at 

where v(A,(/f>,f) is horizontal velocity, h(A,tfi,t) is the thick- 
ness of the fluid layer, and g is gravity. Note that, because of 
the constant-density assumption, the temperature is undefined 
and there is no ther modynamic-energy eq uation. When writ- 
ing these equations, Lan gton & Laughliq ( 1200 7) replaced gh 
with RT, causing the continuity equation (Eq. [TTJ to resem- 
ble an energy equation. This is not a valid procedure, how- 
ever, since the shallow-water layer thickness is a representa- 
tion of mass per area between two bounding material surfaces 
and is thus a distinct quantity from temperature; 11 as men- 
tioned above, the constant-density assumption prevents tem- 
perature from entering the shallow-water system. Their tem- 
perature fields ar e therefore best interp r eted a s shallow-water 
layer thickness. lLangton & Laughlinl ([2007 ) forced Eq. QT| 
using a Newtonian-relaxation scheme that adds mass on the 
dayside and removes it on the nightside, which is intended 
as a simple means to represent dayside radiative heating and 
nightside cooling. Their simulations exhibit development of 

winds reaching ~ lkmsec - ^ 

In a subsequent study, Langton & Laughlin (2008) per- 
formed global, 2D hydrodynamic simulations of the flow on 
hot Jupiters with eccentric orbits. This system of equations, 
unlike the shallow-water set, contains independent continu- 
ity and energy equations; governing variables are v, p, and 
T over the globe. They treated radiation using a two-band 
model, one for stellar irradiation and the other for planetary 
thermal radiation under the assumption that the local emis- 
sion occurs as a blackbody. Their simulations develop strong 
lateral temperature gradients as one side of the planet is flash 
heated during periastron passage; this leads to complex and 
highly turbulent flow fields as the resulting vortex structures 
become dynamically unstable. This natural development of 
turbulence from the hemispheric-scale stellar forcing is an 
interesting result that has not occu rred to date in models of 
hot Jupiters with zero eccentricity ([Showman & Guillot 2002; 
ICoooer & Showmanll200l |2006; Langton & Laughlin) l2007 : 
IDobbs-Dixon & Linll2008l) . and i t may have important impli - 
cations for infrared light curves (lLangton & Laughlinl 2008). 
Nevertheless, 2D dynamics differs in important ways from 3D 
dynamics, and it remains to be demonstrated whether their 
adopted 2D model can reproduce the known flow fields on 

jet structure relative to both their reduced-opacity cases and our results. These 
differences presumably result from the extreme opacity in these simulations. 

11 Equating gh to RT is equivalent to assuming h = RT/g, which is just 
the atmospheric pressure scale height derived for a compressible ideal-gas 
equation of state. This is incompatible with the incompressibility assumption 
on which shallow-water is based; furthermore, a scale height is defined as the 
height between two isobars whereas the shallow-water layer thickness is the 
height between two surfaces of constant potential-density (e.g., density in the 
ocean and entropy in the atmosphere). 
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solar-system planets, which is an important benchmark for 
any model. Exploring 3D circulation models of hot Jupiters 
on eccentric orbits is clearly an interesting avenue for future 
re search. 

ICho et alJ (120031 120081) performed global one-layer sim- 
ulations of the flow on hot Jupiters using the equivalent- 
barotropic equations, which are essentially a vertically aver- 
aged version of the 3D primitive equations. A ma jor differ- 
ence with th e other studies described here is that ICho et alJ 
(20031 120081) did not include any radiative heating/cooling; 
instead, they forced their flows using a combination of small- 
scale turbulence in the initial condition and a hemispheric- 
scale pressure deflection intended to qualitatively represent 
the dynamical effects of the day-night heating gradient. In 
cases when strong initial turbulence was included, their flows 
developed meandering polar vortices, waves, and a high de- 
gree of turbulent mixing. The turbulent initial condition that 
enables these outcomes was intended as a parameterization of 
turbulence generation by atmospheric instabilities. Although 
all planetary atmospheres are turbulent, the appropriate turbu- 
lent length scales and energetic amplitudes remain unknown 
for hot Jupiters, and may well vary from planet to planet. If 
such turbulence is indeed present at the levels explored by 
Cho et al., it should ideally develop naturally in models that 
force the flow (for example with radiative heating/cooling). 
So far, however, radiatively forced investigations of hot 
Jupiters in circular orbits have exhibited relatively steady cir- 
culation patterns that lack the degree of turbulent variabil- 
ity seen in the Cho et al. simulations (Showman & Guilloj 
2002t ICoooer & Showmanll2005t ILaiigton & Laughlinll2007l: 
Dobbs-Dixon & Lira |2008[ and the present stud y), with the 
notable exception of lLangton & Laughlinl |2008)'s 2D study 
of hot Jupiters in highly eccentric orbits. It will be interesting 
to see if this continues to be the case as a wider range of hot 
Jupiters is explored and model resolutions and forcing realism 
improve over time. 

I n contrast to the above dynami cal models. Ilro et al.l ([2005) 
and iBurrows et al.l d2006l [2007b, a) explored models of hot 
Jupiters that used realistic radiative-transfer schemes but 
made ad hoc assum ptions fo r the e ffect of dynamics on 
the thermal budget. Ilro et alJ (120051) assumed atmospheric 
winds that correspond to solid-body rotation; the speed at 
the equator wa s a free parame ter varied from 0.5-2 km sec" 1 . 
IBurrows et al. ( 2006. l2007blah simply removed an arbitrary 
fraction of the stellar irradiation from the dayside and placed 
it on the nightside to mimic the effect of atmospheric winds. 
This is a novel approach that allows the exploration of how 
(for example) day-night temperature differences should de- 
pend on the day-night heat transport. As these authors are 
careful to point out, however, it is not a rigorous treatment of 
dynamics. 

6. THE VALIDITY OF HYDROSTATIC BALANCE 

A scaling analysis demonstrates that local hydrostatic bal- 
ance is approximately valid for the large-scale flow on 
hot Jupiters. It is important to emphasize that the local- 
hydrostatic -balance assumption in the primitive equations de- 
rives from the assumption of large aspect ratio and not from 
any assumption on wind speed. However, the fact that esti- 
mated wind speeds on hot Jupiters are several km sec" 1 , which 
is close to the 3-kmsec" 1 speed of sound in these atmospheres, 
suggests that we consider the validity of the primitive equa- 
tions in the hot- Jupiter context. The full Navier-Stokes verti- 



cal momentum equation can be written 

dw 1 dp 

— + v-Vw = — -^--g + 2uncos(f> (12) 
at p dz 

where w is vertical wind speed, t is time, v is the horizon- 
tal wind velocity, and u is the east-west wind speed. The 
background static hydrostatic balance is irrelevant to atmo- 
spheric circulation and can be removed from the equation. 
Define p = po(z) + p' and p = po(z) + p' where po and po are 
the time-independent basic-state pressure and density and, by 
construction, dpo /dz = -pog- Primed quantities are the devi- 
ations from this basic state caused by dynamics. Substituting 
these expressions into Eq. [12] we can rewrite the equation as 

v+V'Vw = -- (^--p'g) +2uQcoscj> (13) 
at p \ dz J 

where the basic-state hydrostatic balanced has been sub- 
tracted off. The terms in parentheses give only the flow- 
induced contributions to vertical pressure gradient and weight 
(they go to zero in a static atmosphere). 

For local hydrostatic balance to be a reasonable approxi- 
mation, the terms in parenthesis must be much larger than 
the other terms in the equation. The magnitude of dw/dt is 
approximately w/t, where r is a flow evolution timescale, 
and the magnitude of v • Vw is the greater of uw/L and 
w 2 /H where L is the horizontal flow scale and H is the 
vertical flow scale. For global-scale flows, L ~ 10 7 -10 8 m 
and r ~ 10 4 -10 5 sec. The large-s cale flow varies vertically 
over a scale height H ~ 300 km {Showman & Guillot j l2002h 
ICooper & Showmanll200l 120061; IDobbs-Dixon & Lin[[2008l). 
For th e si mulated flow regime inlCooper & S howman (2005, 
120061) and IDobbs-Dixon & Lid {2008J), v ~ k ~ 3 km sec" 1 , 
g ~ 10msec" 2 , fl ~ 2 x 10" 5 sec _1 , and w ~ 10-lOOmsec" 1 . 
With these values, we find that dw/dt < 10~ 2 msec" 2 , uw/L< 
0.03msec" 2 , w 2 /H < 0.03msec" 2 , and Vlu ~ 0.1msec" 2 . In 
comparison, for a hot Jupiter with day-night temperature dif- 
ferences of several hundred K, the flow-induced hydrostatic 
terms p'g/p and p~ l dp' /dz are each ~ 10msec" 2 . 

This analysis implies that, for global-scale hot- Jupiter flows 
at the atmospheric pressures considered in our model (p > 
lmbar), the greatest departure from hydrostaticity results 
from the vertical Coriolis force, which causes a ~ 1% devi- 
ation from hydrostatic balance. The acceleration terms on the 
left side of Eq. Qj] (which are necessary for vertically prop- 
agating sound waves) cause a ~ 0.3% deviation from hydro- 
static balance. Hydrostatic balance is thus a reasonable ap- 
proximation for the large-scale flow. Nonhydrostatic effects 
of course become important at small scales, and it is conceiv- 
able that these effects interact with the large-scale flow in non- 
trivial ways. For hot Jupiters, the acceleration terms on the 
left side of Eq. [13] only become important for structures with 
vertical or horizontal scales less than ~ 30 km and ~ 500- 
1000km, respectively. In a numerical model that solved the 
full Navier-Stokes equations, the grid resolution would have 
to be substantially finer than these values for the nonhydro- 
static behavior to be accurately represented. A full Navier- 
Stokes solution with a coarse resolution would effectively be 
resolving just the global-scale hydrostatic component of the 
flow. 

Despite the above, some hot Jupi ters may be losing mass 
from the top of their atmospheres (Vidal-Madiar et aill2003l 
2004), which suggests that hydrostatic equilibrium should 
break down at extremely low pressures (above the top of our 
model) where this outflow occurs. 
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We wish here to clarify the issue of vertical velocities. 
iDobbs-Dixon & LirJ (120081) suggested that vertical motions 
are faster in the full Navier-Stokes system than in the primitive 
equations because of the different vertical momentum equa- 
tions in these systems. Certainly, velocities in a convection 
zone require use of a momentum equation that includes ver- 
tical accelerations. However, in statically stable atmospheres, 
the primary control over the vertical velocity in a quasi-steady 
overturning circulation (relevant to the regime simulated here) 
is not the vertical momentum equation but the thermodynamic 
energy equation (Eq. |4j. The reason is that, in this situation, 
the radiative heating/cooling rate controls the speed of vertical 
motion. Because entropy increases with height in a statically 
stable atmosphere, adiabatic expansion/contraction in ascend- 
ing (descending) air causes temperature at a given height to 
decrease (increase) over time. In the absence of radiation, 
such steady flow patterns are unsustainable because they in- 
duce density perturbations that resist the motion (i.e. ascend- 
ing air becomes denser and descending air becomes less dense 
than the surroundings at that altitude). Thus, steady vertical 
motion in a stable atmosphere can only occur as fast as radi- 
ation can remove the temperature perturbations caused by the 
adiabatic ascent/descent. This provides a fundamental con- 
straint on vertical motion that applies equally to the primitive 
and the Navier-Stokes equations. 

To quantify this idea, consider Eq. |U which can be written 



dT H 2 N 2 q 

— +v vr-w = — 

at Rp c p 



(14) 



where H is the vertical scale height and N is the Brunt- Vaisala 
frequency (the oscillation frequency for a vertically displaced 
air parcel in a statically stable atmosphere). An upper limit on 
the attainable vertical velocity in a quasi-steady circulation in 
a statically stable atmosphere results from equating the right 
side to the third term on the left side. One then obtains a peak 
vertical velocity (measured in Pa sec" 1 ) 



q Rp 

c p H 2 N 2 



(15) 



which can be approximately expressed as a peak vertical ve 
locity in msec -1 

q R 



c p HN 2 



(16) 



The equation of course couples to the rest of the dynam- 
ics via Af 2 and q/c p . For hot Jupiter photospheres, where 
q/cp ~ 10 -2 Ksec _1 , H ~ 300km, and N ~ 0.001-0.01 sec" 1 
depending on temperature, gravity, and thermal gradient, we 
then e xpect w ~ 10-lOOmsec -1 or less (Sh owman & Guillol 
2002). Importantly, Eq. [T6l makes no assumption about the 
form of the vertical momentum equation (hydrostatic versus 
Navier-Stokes); only quasi-steady flow and a statically stable 
atmosphere were invoked. Within the context of the static sta- 
bilities and heating rates assumed above, these vertical veloc- 
ities should therefore apply to statically stable atmospheres 
regardless of whether they are modeled with the primitive 
or Navier-Stokes equations. On the other hand, because of 
the intense nightside cooling, some hot Jupiters may develop 
convection zones near the photosphere on the nightside. The 
primitive equations cannot capture the small-scale convective 
vertical velocities in this situation. Nevertheless, our simu- 
lations of HD209458b and HD 189733b appear not to predict 
such detached convection regions. 



7. CONCLUSIONS 

We presented three-dimensional numerical simulations of 
the atmospheric circulation on HD209458b and HD 189733b; 
calculated infrared spectra and light curves from the simu- 
lated temperature patterns; and compared these observables 
with available measurements. The simulations were forced 
with a simplified Newtonian heating/cooling scheme. Within 
the context of this approach, we have improved on earlier 
work by calculating for use in this scheme the true radiative- 
equilibrium temperatures as a function of longitude, latitude, 
and pressure; likewise, we calculated the radiative time con- 
stants as a function of both temperature and pressure over the 
3D grid. These radiative ti me constants are generally consis- 
tent with earlier estimates dlro et al.l 12005b but span a much 
greater range of pressure and especially temperature. 

In our simulations, HD209458b and HD 189733b develop 
similar circulation patterns, although the latter has a global- 
average temperature ~ 200-300 K cooler than the former (a 
result of the lower stellar flux). C onsistent with earlier work 
(Cooper & Showman 2005, 2006), our simulations show that, 
at low pressure (< lOmbar), the circulation transports air 
from dayside to nightside, both along the equator and over 
the poles. At high pressure (> lOOmbar), a banded structure 
emerges with a broad, fast (3-4 km sec -1 ) eastward equatorial 
jet flanked by weaker westward flow at high latitudes. The 
day-night temperature difference varies strongly with height 
and reaches 500-1000K above the photosphere. Near the 
photosphere, the dynamics distorts the temperature pattern 
in a complex manner; consistent with our earlier work, this 
generally includes an eastward displacement of the hottest 
regions from the substellar point. Importantly, dynamics 
also pushes the vertical temperature profile Tip) far from 
radiative-equilibrium. The temperature decreases strongly 
with altitude on the nightside but becomes quasi-isothermal, 
or even exhibits an inversion layer, on the dayside. 

Our calculated spectra, calculated assuming a cloud-free at- 
mosphere with solar metallicity (appropriately modified for 
rainout) show that the nightside should exhibit deep absorp- 
tion bands. On the dayside, however, the deviations of the 
spectrum from a blackbody become modest due to the quasi- 
isothermal structure. Nevertheless, in some simulations this 
dayside structure contains weak emission features, which re- 
sult from the existence of a thermal inversion in these simula- 
tions. 

Our light curves calculated in Spitzer bandpasses show 
that, for HD 189733b, we correctly explain the phase offset 
of the flux peak that occurs befor e the secondary eclips e mea- 
sured in the 8-^m light curve of lKnutsonetail ([2007b), but 
we fail to explain the flux minimum that occurs after transit 
and we overpredict the total flux variation. For HD209458b, 
we match the Spitzer IRAC 3.6 and 8-^m secondary-eclipse 
depths, and marginally fit the latter of two 24 fim MIPS mea- 
surements, but we underpredict the 4.5 and 5.8 /im IRAC 
secondary-eclipse depths. This probably results from the lack 
of a strong stratosphere that seems to be implied by these ob- 
servations. Future simulations that include a realistic repre- 
sentation of radiative transfer, including the possibility of TiO 
and VO opacity for the case of HD209458b, may be needed 
to better fit these diverse observations. 

Finally, we presented a scaling analysis suggesting that the 
large-scale (> 10 7 m) flow on the hot Jupiters explored here 
should be close to local hydrostatic balance; on these large 
scales, deviations from hydrostaticity are small. This re- 
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suit lends support to the primitive-equation approach adopted 
here. At the present time, the primary source of uncertainty is 
not the dynamical scheme but the representation of the forcing 
that drives the flow. 

So far, all published studies of the circulation on hot 
Jupiters have forced the flow using relatively severe approx- 
imations to the radiative transfer or excluded heating/cooling 
entirely. Inclusion of realistic radiative transfer is an impor- 
tant goal for future work. Other areas for improvement in- 
clude considering cloud/haze formation, alternate composi- 
tions (supersolar metallicity and the possibility of TiO and 
VO), and exploring a wider range of planetary parameters to 



capture the diversity of the growing number of known hot 
Jupiters. 
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Radiative Time Constants for HD 209458b 
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NOTE. — Column headers are temperatures in Kelvin. Values are time constants in seconds. Note that no values were calculated in 
the lower left corner here and in Table 2 because the radiative-equilibrium p-T profiles used to calculate r ra[ ] (Fig.fl} do not access these 
pressure/temperature combinations. 



TABLE 2 

Radiative Time Constants for HD 189733b 
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Note. — 



Column headers are temperatures in Kelvin. Values are time constants in seconds. 



